Vegetation Change¶

This notebook uses changes in the Normalized Difference Vegetation Index (NDVI) to identify vegetation change. The algorithm identifies a "baseline" and "analysis" time period and then compares the spectral index in each of those time periods. Significant changes in NDVI (vegetation greenness) are coincident with land change, as long as the comparisons are done between similar time periods (seasons or years). Users of this algorithm should not accept the accuracy of the results but should conduct ground validation testing to assess accuracy. It is expected that this algorithm can be used to identify clusters of pixels that have experienced change and allow targeted investigation of those areas by local or regional governments. In some cases the impacts may be negative (deforestation, mining, burning, drought) or positive (regrowth, improved soil moisture).

It should also be noted that the selection of the baseline and analysis time period is critical. First, the two time periods should be similar (season, year) so that the vegetation state can be compared in similar weather conditions. Second, the time periods should be sufficiently clear (non-cloudy) data. If the baseline or analysis mosaic (composite of images) is contaminated with clouds, it will impact the results.

Connect to the Data Cube ▴¶

In [1]:
import datacube
dc = datacube.Datacube(app='Vegetation_Change')

import sys, os
os.environ['USE_PYGEOS'] = '0'

from datacube.utils import masking
from dea_tools.plotting import rgb, display_map
from odc.algo import to_f32

### EASI tools
sys.path.append(os.path.expanduser('../scripts'))
from ceos_utils.data_cube_utilities.clean_mask import landsat_clean_mask_invalid, landsat_qa_clean_mask
from easi_tools import EasiDefaults
from easi_tools import notebook_utils
easi = EasiDefaults() # Get the default parameters for this system

Successfully found configuration for deployment "eail"
In [2]:
cluster, client = notebook_utils.initialize_dask(use_gateway=False)
display(cluster if cluster else client)
print(notebook_utils.localcluster_dashboard(client, server=easi.hub))

LocalCluster

64a0b9b9

Dashboard: http://127.0.0.1:8787/status Workers: 4
Total threads: 8 Total memory: 29.00 GiB
Status: running Using processes: True

Scheduler Info

Scheduler

Scheduler-1d3a4938-5207-40b6-aa9c-3b0bbb6c36c3

Comm: tcp://127.0.0.1:43369 Workers: 4
Dashboard: http://127.0.0.1:8787/status Total threads: 8
Started: Just now Total memory: 29.00 GiB

Workers

Worker: 0
Comm: tcp://127.0.0.1:42903 Total threads: 2
Dashboard: http://127.0.0.1:39163/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:40295
Local directory: /tmp/dask-worker-space/worker-twl206p5

Worker: 1
Comm: tcp://127.0.0.1:37309 Total threads: 2
Dashboard: http://127.0.0.1:42975/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:42435
Local directory: /tmp/dask-worker-space/worker-7cansc5h

Worker: 2
Comm: tcp://127.0.0.1:38089 Total threads: 2
Dashboard: http://127.0.0.1:38395/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:44559
Local directory: /tmp/dask-worker-space/worker-pz6usymu

Worker: 3
Comm: tcp://127.0.0.1:37713 Total threads: 2
Dashboard: http://127.0.0.1:37015/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:46389
Local directory: /tmp/dask-worker-space/worker-damtr0hx 
https://hub.eail.easi-eo.solutions/user/jhodge/proxy/8787/status
In [3]:
from datacube.utils.aws import configure_s3_access
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client)
Out[3]:
<botocore.credentials.Credentials at 0x7f048557d780>

Choose Platforms and Products ▴¶

In [4]:
# Select a Product
product = "landsat8_c2l2_sr"

Define the Extents of the Analysis ▴¶

In [5]:
# Select an analysis region (Latitude-Longitude) 

# Mining Region near Obuasi, Ghana 
# Baseline = 12/23/13, Analysis = 1/6/19 
# NOTE: These are clear scenes at similar times of the year. 
# Use the Cloud Statistics notebook to evaluate cloud cover.

latitude = (-23.5,-23.6)
longitude = (-59.45,-59.35)
In [6]:
## The code below renders a map that can be used to orient yourself with the region.
display_map(longitude,latitude)

/env/lib/python3.10/site-packages/dea_tools/plotting.py:313: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
  all_longitude, all_latitude = transform(Proj(crs), Proj('EPSG:4326'), all_x,
Out[6]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Define Analysis Parameters ▴¶

In [7]:
from datetime import datetime

# Select the start and end periods for your analysis products (Year,Month,Day).
# These time windows will be used to make a median mosaic, so typically pick a year length (or more)
# or select a small window surrounding a clear single date (use Cloud Statistics notebook).
# Also, be sure to evaluate the RGB mosaics (below) to affirm they are not full of clouds.

# Select the baseline time period (start and end)
baseline_time_period = (datetime(2013,11,12), datetime(2013,11,13))

# Select the analysis time period (start and end)
analysis_time_period = (datetime(2022,12,23), datetime(2022,12,24))

Load and Clean Data from the Data Cube ▴¶

In [8]:
output_crs = 'epsg:6933'
resolution = (-30,30)
scale = 0.0000275
offset = 0.2
In [9]:
baseline_ds = dc.load(latitude=latitude,longitude=longitude,product=product,
                      measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'],
                      time=baseline_time_period,output_crs=output_crs,resolution=resolution,group_by='solar_day',dask_chunks={'time':1})
In [10]:
analysis_ds = dc.load(latitude=latitude,longitude=longitude,product=product,
                      measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'],
                      time=analysis_time_period,output_crs=output_crs,resolution=resolution,group_by='solar_day',dask_chunks={'time':1})
In [11]:
cloud_mask_baseline = masking.make_mask(baseline_ds['pixel_qa'], clear='clear', cloud='not_high_confidence', cloud_shadow='not_high_confidence')
baseline_ds = baseline_ds.where(cloud_mask_baseline)

cloud_mask_analysis = masking.make_mask(analysis_ds['pixel_qa'], cloud='not_high_confidence', cloud_shadow='not_high_confidence')
analysis_ds = analysis_ds.where(cloud_mask_analysis)

Create a Median Mosaic for the Baseline and Analysis Time Periods¶

In [12]:
from ceos_utils.data_cube_utilities.dc_mosaic import create_median_mosaic
baseline_composite = to_f32(create_median_mosaic(baseline_ds, cloud_mask_baseline), scale=scale, offset=offset)
analysis_composite = to_f32(create_median_mosaic(analysis_ds, cloud_mask_analysis), scale=scale, offset=offset)

Calculate NDVI Anomaly Product ▴¶

In [13]:
def NDVI(dataset):
    return (dataset.nir - dataset.red)/(dataset.nir + dataset.red)
In [14]:
parameter_baseline_composite = NDVI(baseline_composite)
parameter_analysis_composite = NDVI(analysis_composite)
parameter_anomaly = parameter_analysis_composite - parameter_baseline_composite

Show Baseline RGB, Analysis RGB and Anomaly Products¶

In [15]:
import matplotlib.pyplot as plt
from ceos_utils.data_cube_utilities.dc_rgb import rgb
from matplotlib.cm import RdYlGn
RdYlGn.set_bad('black',1.)
In [16]:
# Define the significant anomaly range for Plot #4
# The typical loss range is <0. Choose < -0.1 or < -0.2.
# The typical gain range is >0. Choose > 0.1 or > 0.2. 

loss_range  = parameter_anomaly < -0.05
gain_range  = parameter_anomaly > 0.1
In [17]:
import xarray as xr
import numpy as np

fig, ax = plt.subplots(2, 2, figsize=(12,12))

for sub_ax in ax.flatten():
    sub_ax.set_facecolor('black')

baseline_rgb = baseline_composite[['red', 'green', 'blue']].to_array().compute()
analysis_rgb = analysis_composite[['red', 'green', 'blue']].to_array().compute()

# Use the middle values of the data (2% to 98%) to brighten the image
lw_qtl, up_qtl = 0.02, 0.98
rgb_vmin = min(baseline_rgb.quantile(lw_qtl).values,analysis_rgb.quantile(lw_qtl).values)
rgb_vmax = max(baseline_rgb.quantile(up_qtl).values,analysis_rgb.quantile(up_qtl).values)

# Plot the resulting 4 products ... Baseline RGB, Analysis RGB, Total Anomaly, Anomaly Threshold
# NOTE: Clouds in either the baseline or analysis images will be removed from the anomaly product

## Plot #1 = Baseline RGB (upper left)
axes_image = baseline_rgb.plot.imshow(ax=ax[0,0], vmin=rgb_vmin, vmax=rgb_vmax)

## Plot #2 = Analysis RGB (upper right)
analysis_rgb.plot.imshow(ax=ax[0,1], vmin=rgb_vmin, vmax=rgb_vmax)

## Plot #3 = Total Anomaly (lower left)
parameter_anomaly.plot(ax=ax[1,0], vmin=-0.2, vmax=0.2, cmap = RdYlGn, add_colorbar=False)

## Plot #4 = Anomaly Threshold (lower right)

# Analysis composite grayscale background
plt4_bkg_band = 'swir1' # The band to use as the background image.
plt4_rgb = np.repeat(analysis_composite[plt4_bkg_band].where(cloud_mask_baseline.sum('time').astype(np.bool))\
                     .values[:,:,np.newaxis],3,axis=2)

# Selected a range of SWIR1 values (0.001 to 0.600) to lighten image background 
# Users may also try values of 0.02 and 0.98
min_bkg = np.nanquantile(analysis_composite[plt4_bkg_band].values, 0.001)
max_bkg = np.nanquantile(analysis_composite[plt4_bkg_band].values, 0.95)
plt4_rgb = np.interp(plt4_rgb, (min_bkg, max_bkg), [0,1])

# Significant anomaly color overlays
color_green = np.array([0,1,0]) # green
color_red   = np.array([1,0,0]) # red
plt4_rgb[loss_range] = color_red
plt4_rgb[gain_range] = color_green

# Plot
plt4_coords = dict(analysis_composite.coords)
rgb_coord_arr = np.array(['red', 'green', 'blue'])
rgb_coord_da = xr.DataArray(rgb_coord_arr,name='rgb',dims=['rgb'],coords={'rgb': rgb_coord_arr})
plt4_coords.update({'rgb': rgb_coord_da})
plt4_rgb_da = xr.DataArray(plt4_rgb, coords=plt4_coords,dims=list(analysis_composite.dims) + ['rgb'])
plt4_rgb_da.plot.imshow(ax=ax[1,1])

# Titles for all plots
ax[0,0].set_title('Baseline Composite'), ax[0,0].xaxis.set_visible(False), ax[0,0].yaxis.set_visible(False)
ax[0,1].set_title('Analysis Composite'), ax[0,1].xaxis.set_visible(False), ax[0,1].yaxis.set_visible(False)
ax[1,0].set_title('Vegetation Anomalies: Red=Loss, Green=Gain'), ax[1,0].xaxis.set_visible(False), ax[1,0].yaxis.set_visible(False)
ax[1,1].set_title('Locations of Significant Anomalies: Red=Loss, Green=Gain'), ax[1,1].xaxis.set_visible(False), ax[1,1].yaxis.set_visible(False)

plt.tight_layout()
plt.show()

/tmp/ipykernel_510/2898797902.py:33: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  plt4_rgb = np.repeat(analysis_composite[plt4_bkg_band].where(cloud_mask_baseline.sum('time').astype(np.bool))\
In [ ]: